
    CalMarRatio = function(pars, data, cnt, eta.coef, rep,seed){
        
        ############# part 1: get density of L0 #########################
        
        L0 = cbind(rep(1,nrow(data)),data[,1:pL0])
        l0 = L0 %*% 2^(1:(pL0+1))
        L0 = unique(L0,margin = 1)
        l0.index = match(l0, unique(l0))
        Nb = nrow(L0)
        density.l0 = as.numeric(by(cnt,l0.index,sum)) / sum(cnt)
        
        
        ############# part 2: get E[Y|\bar{A}_K=1, \bar{L}_K, L_0] ##########
        
        theta.coef = matrix(pars$theta.coef[c(1,3:(pL0+2), 2:(pL0+2))],,2)
        colnames(theta.coef) = c("l0","l1")
        theta = exp(L0 %*% theta.coef)
        
        phi.coef   = matrix(pars$phi.coef[c(1,4:(pL0+3),2,4:(pL0+3),
                                            3:(pL0+3))],,3)
        colnames(phi.coef)   = c("a0","a1","a.") 
        phi   = exp(L0 %*% phi.coef)
        eta   = expit(L0 %*% eta.coef)
        gop   = exp(L0 %*% pars$gop.coef)
        
        ### Step 1: get ratio between two "adjacent" probabilities
        ratio = t(sapply(1:Nb, function(i) GetRatio(theta[i,], phi[i,],eta[i,])))
        dim(ratio) = c(Nb, 5, 2)
        dimnames(ratio) = list(1:Nb, c("ak","lk","aK","lK","l0"),c("0","1"))
        
        ### Step 2: get the maximal ratio against the "baseline p000" 
        #   for each observed unit (only depends on L0)
        k.max = GetMaxKGen(ratio, K, Nb)
        k.max[k.max == Inf] = .Machine$double.xmax
        check = which(k.max < 1e10)
        
        ### Step 3: get the largest p (only depends on L0)
        p.max = rep(1 - 0.1^10, Nb)
        if(length(check)>0)     
            p.max[check] = PMaxEstRough(ratio[check,,], K, k.max[check],
                                        gop[check],rep = rep, seed = seed)
        
        ### Step 4: get the ratio against the largest p for each observed unit
        L = t(sapply(0:(2^(K+1)-1), function(x) 
            as.integer(intToBits(x))[1:(K+1)]));   colnames(L) = 0:K
        A1 = matrix(1, nrow(L), K+1);   colnames(A1) = 0:K
        A0 = matrix(0, nrow(L), K+1);   colnames(A0) = 0:K
        
        k.a1 = GetKGen2(A1, L, ratio, k.max, K)
        k.a1[k.max == .Machine$double.xmax] = 0
        p.y.a1 = p.max * k.a1
        k.a0 = GetKGen2(A0, L, ratio, k.max, K)
        k.a0[k.max == .Machine$double.xmax] = 0
        p.y.a0 = p.max * k.a0
        
        
        
        ############## part 3: estimate P(L_k+1|\bar{A}_k=a, \bar{L}_k, L_0) ##
        
        GetProb  = function(p, x0)    (p^x0) * ((1-p)^(1-x0))
        GetProb2 = function(p0, p1, xk, xk_1){
            (p0^((1-xk) * xk_1)) * ((1-p0)^((1-xk) * (1-xk_1))) *
            (p1^(xk * xk_1))     * ((1-p1)^(xk * (1-xk_1)))
        } 
        
        p.L.a1 = p.L.a0 = outer(eta[,"."],L[,"0"], GetProb)
        for(k in 1:K){
            for(i in 1:nrow(p.L.a1)){
                p.L.a1[i,] = p.L.a1[i,] * GetProb2(eta[i,"10"], eta[i,"11"],
                                                   L[,k], L[,k+1])  
                p.L.a0[i,] = p.L.a0[i,] * GetProb2(eta[i,"00"], eta[i,"01"],
                                                   L[,k], L[,k+1]) 
            }
        }
        
        EY.a1 = sum(rowSums(p.y.a1 * p.L.a1) * density.l0)
        EY.a0 = sum(rowSums(p.y.a0 * p.L.a0) * density.l0)
        
        return(c(EY.a1, EY.a0, EY.a1/EY.a0))
        
    }
    
    
    
    
    
    
    
    
    
    